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Electrolytes Between Dielectric Charged Surfaces: Simulations and Theory 
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We present a simulation method to study electrolyte solutions in a dielectric slab geometry using a modified 
3D Ewald summation. The method is fast and easy to implement, allowing us to rapidly resum an infinite 
series of image charges. In the weak coupling limit, we also develop a mean-field theory which allows us to 
predict the ionic distribution between the dielectric charged plates. The agreement between both approaches, 
theoretical and simulational, is very good, validating both methods. Examples of ionic density profiles in 
the strong electrostatic coupling limit are also presented. Finally, we explore the confinement of charge 
asymmetric electrolytes between neutral surfaces. 


I. INTRODUCTION 

Complex charged liquids present a variety of inter¬ 
esting phenomena such as like-charge attractioni*^ and 
charge inversion^^— . These phenomena arise as a result 
of electrostatic correlations of counterions in the dou¬ 
ble layer—. To explore the ionic correlations between the 
overlapping double layers one must have accurate and 
reliable methods to study electrolytes in a slab geome¬ 
try. Similar difficulties are encountered when one wants 
to understand the thermodynamics of ionic liquid based 
supercapacitorsi. 

Electrolytes in a slab-like geometry have been exten¬ 
sively studied in the past— . However, the dielectric 
contrast between the solvent and the confining surfaces is 
usually not taken into account. This dielectric mismatch 
leads to polarization of interfaces which dramatically in¬ 
creases the mathematical complexity of the problem. For 
a single interface, the dielectric heterogeneity and the re¬ 
sulting induced surface charge can lead to a repulsion of 
electrolyte from the charged dielectric surfac o^^d^ . Sim¬ 
ilar behavior has been observed near charged colloidal 
particlesi^^— . 

Most of the theoretical and simulational works involv¬ 
ing confined electrolytes neglect the dielectic contrast 
which results in an infinite series of image charges. While 
this significantly simplifies the calculations, it also fails 
to account for some of the fundamental physics of the 
overlapping double layers. Recently Wang and Wang22 
presented a mean-field theoretical discussion of confined 
electrolytes between charged and neutral plates. In the 
same year, Zwanikken and de la Crua^ developed a liq¬ 
uid state theory which predicted that neutral confining 
polarizable interfaces can attract each other inside an 
electrolyte solution. Similar result can be found in other 
works2i“— . Samaj and Trizao^l and Jho et al^ de¬ 
veloped theories to study the distribution of confined 
counterions between charged plates, in a salt free sys¬ 
tem in a strong-coupling limit. Some studies focused on 
the specific case in which the dielectric constant of the 
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surrounding medium is much lower than of water— 
Jho et al^ developed a simulation method for con¬ 
fined counterions based on the electrostatic layer cor¬ 
rection (EEC) method^. Also, Tyagi et con¬ 

structed the ICMMM2D method which is the generaliza¬ 
tion of the MMM2D algorithni^, previously developed 
to study homogeneous dielectric slab systems. Although 
these methods account for the surface polarization, they 
require a calculation of a sum of terms for the electro¬ 
static potential that make simulations quite slow. Simi¬ 
lar difficulties are encountered with other simulation ap- 
proaches^^. Boundary element methods (BEM) consider 
the minimization of functionals and can be applied to 
systems with general geometries. Some BEM methods 
consider local polarization charge densities as dynamic 
variables^i^^ others, attempt to explicitly calculate the 
bound charge^SiiS. Even though BEM methods are ex¬ 
pensive computationally, they have been extensively used 
to study general soft matter problems^^— . 

In this paper we present a simulation method based on 
3D Ewald summation with a modified correction for the 
slab geometry^. The method is simple to implement. It 
does not require summations of a slowly convergent infi¬ 
nite series of images during the simulation, and is com¬ 
parable in time with a regular 3D Ewald method. The 
paper is organized as follow. In Section |B1 we show how 
to construct the electrostatic energy of the system. In 
Section uni we study confined electrolytes between po¬ 
larizable charged surfaces. In Section IIYI we present a 
mean-field theory for 1:1 electrolytes and compare it with 
simulations. In Section |Vl we present the general results, 
and in the Section ED the conclusions. 


II. ELECTROSTATIC ENERGY 

To perform simulations, we use a rectangular simu¬ 
lation box of sides Ly and L^. The box contains 
N ions of charges qj = ajQ, where aj is the valence 
of the ion and Q is the proton charge, confined in the 
region —Lxj^ < x < —Lyl2 < y < Ly/2 and 

—L/2 < z < Ll2. We set Ly = and = 5L. The 

uniform dielectric constants are: inside the slab con¬ 
taining electrolyte, and to outside. The dielectric con¬ 
trast results in an infinite set of “images of images” which 
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must be resumed to obtain the total electrostatic energy. 
We define as a number of images of an ion at each 
interface. To calculate the exact electrostatic energy, Ni 
should be infinite. This, however, is not practical in a 
simulation. Instead, we explore the convergence of sim¬ 
ulations as the number of images W is increased. For 
example, if W = 2 , we consider one image charge at each 
dielectric interface and the image of image, producing 4 
image charges for each ion, see Fig. [TJ The electrostatic 
potential at the position r (in the region with e^,), cre¬ 
ated by all ions (excluding ion i), their image charges 
(including the image of ion i), and the periodic replicas 
is 



where pj{s) = qjSis - - rep) and pf^(s) = 7 "^qjS(s - 

rf^—Tep) are the charge densities of ions and their repli¬ 
cas and of dielectric images and their replicas. 

The replication vector is defined as rep = L^rixi: -I- 
Lynyy + LzUzZ andr^ = xjx + yjy+[{-!)""zj ±mL]z. 
The vectors n = {nx,ny,nz), where Ux, Uy and Uz 
are positive and negative integers, represent the infinite 
replicas of the main cell. The constant 7 is defined as 
7 = (ctu — eo)/(eu) + £ 0 ) and the prime on the summa¬ 
tion signifies that j ^ i, when n = (0,0,0). The total 
electrostatic energy of the system is 
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( 2 ) 


i=l 


The energy above is very difficult to calculate because 
of a slow conditional convergence of the series in Eq. ©■ 
To speed up the convergence, we use the Ewald method 
in which the ionic charge is partially screened by placing 
a Gaussian-distributed charge of opposite sign on top of 
each ion^. We then add and subtract an opposite Gaus¬ 
sian charge at the position of each ion and its images, 
Pj{s) and p^^(s), respectively. The potential, Eq. ([T]), 
then becomes 
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FIG. 1. Representation of a charge between the dielectric 
surfaces. Only the first and second order images are shown, 
Ni = 2 . 
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Pf(s) = exp (-K^ls-rj - repp) , 


( 5 ) 
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( 7 ) 


pfm (s) = 7™9i(Ke/^) exp {-kI\s - rf^ - repp) , 

( 8 ) 

and Ke is a damping parameter. Note that we have sub¬ 
tracted the self potential, Eq. from Eq. @, in order 


Mr) = Mir) + Mir) - (t>MHr) 
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to remove the prime over the summation in the long- 
range (L) part of the potential, Eq. ([S]). The electro¬ 
static potential produced by the Gaussian charges can 
be calculated using the Poisson equation, yielding 


OO N 

n j^l 

Ni OO N 
m—1 n j—1 


eri{Ke\r-rj -rep\) 
ew\r-rj -Vep] 

eTf{Ke\r-r+^-rep\) 
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eTi{Ke\r -r^^-rep\) 
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(9) 


where erf(a;) is the error function. The short-range part of 
the potential (5), Eq. (jd]), can then be obtained in terms 
of the complementary error function, erfc(a:) = 1 —erf(x), 
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This potential decays exponentially fast, with the decay 
length controlled by the damping parameter which we set 
to Ke = 'i/Rc: where Rc = is the distance cutoff. It 
is then necessary to consider only the term n = (0, 0,0), 
with the usual periodic boundary condition. Further¬ 
more, for sufficiently large values of k its is sufficient to 
include only a few images-of-images. The damping pa¬ 
rameter, however, can not be too high since its value 
controls the number of fc-vectors that will have to be 
used to calculate the long-range potential. For systems 
studied in this paper, we find that Ni = 2, in the short 
range potential, is sufhcient. Depending on the separa¬ 
tion between the plates, more images may be necessary. 
Prior to accumulation of data we, therefore, check for 
convergence by varying the value of Ni . The short-range 
potential then becomes 
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The self-potential, Eq. ®, reduces to 
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We next calculate the long-range part of the potential, 
Eq. ([9]). This is most easily obtained using the Fourier 
representation, (j)^{k) = J dPr exp{—ik-r)(f)^{r), since 
in the reciprocal space all sums, once again, con¬ 
verge very rapidly. The Fourier transform p^{k) = 
f (Pr exp {—ik ■ r)p^(r), of the Gaussian charge density. 
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where k = {2TTnx/Lxy,2TTny/Lxy,2TTnz/Lz) and V = 
LxLyLz- Using the Poisson equation, |A:p(^’"(fc) = 
{4:TT/ew)p^{k), we can evaluate the Fourier transform of 
the potential, 
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The corresponding real-space electrostatic potential is 
calculated using the inverse Fourier transform, (t)^{r) = 

<^’'(A:)exp(tA:-r), 

,r. . dvr . I^P , , 
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The long-range contribution to the total electrostatic 
energy is then given hy Ul = (1/2) gi'/^(^i)i where 

4>^{r) is obtained from Eq. ((T5)) . It is convenient to 
rewrite this in terms of functions: 
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N 

B{k) =-'^qism{k-ri) , (18) 

2=1 


Ni N 

C{k) =J2Y1 [cos {k ■ r+^) + cos {k ■ r~^)] (19) 

m—1 i—1 

and 


Ni N 


= -J2J2 [®™ ■ ^m) + sin (fc • r,^)] . 

m=l 2=1 

( 20 ) 

The electrostatic energy then becomes, 

TT 


.V\k\^ 


[A{kf + B{kf + A{k)C{k) + B{k)D{k)] . (21) 


The terms in Eqs. (HU) and (l20l) . are multiplied by the 
7 ™ parameter, leading to a converging sum for realis¬ 
tic 7 < 1 parameter values. However, we do not know 
a priori a minimum number of images necessary to ob¬ 
tain an accurate result for the long-range potential. For 
example, we find that for 7 si 0.9, we need Ni = 50 to ob¬ 
tain a good convergence. However, such a large number 
of images makes simulations extremely slow. We note, 
however, that Eqs. (HU and dsni) can be rewritten as, 


N 

c(k) = qi [ci(fe) cos (kxXi + kyiji) cos {kzZi) + 

i=l 

C 2 {k) sin {kxXi + kyPi) sin {k:,Zi)l22) 


and 
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D{k) = - '^qi [di(fc) sin {k^Xi + kyPi - k^zi) + 

i=l 

d 2 {k) sin {kxXi + kyPi -I- k^zt)] .(23) 

where. 


Ni 

ci(fc) = 2 ^ 7 ™ cos {mkzL) , 

m—1 


Ni 

C 2 {k) = 2 ^ cos (mkzL) , 

772=1 


Ni 

di{k) = 2 ^ 7 ™° cos (wofezL) , (24) 

7?2o = l 

Ni 

d 2 {k) = 2 ^ 7 ™® cos (rriefezL) , (25) 

m^=2 


and m are integers: rrio are odd and vie are even. The 
parameters ci(fc), C 2 (fc), di{k) and d 2 {k) can be obtained 
once (up to any desired accuracy) at the beginning of 
the simulation, since they do not depend on the ionic 
positions. The functions, A{k), B{k), C{k) and D{k), 
can now be easily updated for each new configuration in 
a Monte Carlo (MC) simulation. 

The electrostatic energy resulting from the short-range 
part of the potential is Us = ( 1 / 2 ) where 

4>i{r) is given by the Eq. (HU, and the self-energy 
contribution is Useif = ( 1 / 2 ) Sili In the 

limit a; —> 0 , the erf(a;) function vanishes as ( 2 /-y^)a: 
and the self-energy contribution reduces to Useif = 
(Ke/eu,i/^) 9 ?- The total electrostatic interaction 

energy of the ions is given by the expressions above, plus 
the correction needed to account for the slab geometry^. 

Yeh and Berkowita^ found that a regular 3D Ewald 
summation method with an energy correction which ac¬ 
counts for the anisotropic summation of a conditionally 
convergent series in a slab-like geometry can reproduce 
the same results as the 2D Ewald method, with a signifi¬ 
cant gain in performance. For more details on the Ewald 
summation method, applied to different geometries, an 
interested reader can consult Refs 4^“—. Following Yeh 
and Berkowitz and Smitb^ and taking into account the 
dielectric discontinuity and the induced image charges, 
we find the correction for the slab geometry to be 
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where = (—l)™Zj±mL. Using the electroneutrality, 
this expression can be rewritten as 


Uear — 
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where Mz = Ih® magnetization in the z di¬ 

rection. Again, the constant between the brackets can 
be evaluated once at the beginning of the simulation, so 
that Mz can be easily updated in the simulation process. 
The total energy of the system is then 


U =Us + UL + Ueor + Uself ■ (28) 


III. MOIMTE CARLO SIMULATIONS 

We now study an electrolyte solution confined be¬ 
tween two negatively charged dielectric surfaces. The 
two charged plates contain 256 point charges each, uni¬ 
formly distributed on the surface on a square lattice. The 
magnitude of point charges is adjusted to obtain the de¬ 
sired surface charge density. The surfaces are located 
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at z = —Ll2 and z = +L/2. Counterions are modeled 
as hard spheres with the charge q located at the center. 
Besides the counterions, salt ions can also be present in 
the system, all with the same ionic radius, 2 A. The sol¬ 
vent is modeled as an uniform dielectric medium. The 
Bjerrum length As = jtw of the system is 7.2 A, cor¬ 
responding to water at room temperature. The MC sim¬ 
ulations are performed using the Metropolis algorithm. 
The method developed in Section [H] is used to obtain the 
electrostatic energy. Care must be taken in the calcula¬ 
tion of the electrostatic energy of wall particles. For these 
particles the self image electrostatic interaction diverges, 
leading to an infinite constant which must be renormal¬ 
ized. We use 1 X 10® attempted particle moves to equili¬ 
brate and 100 moves per particle to create a new state. 
After 40000 states, we calculate the average ionic density 
profiles. 


IV. MODIFIED POISSON-BOLTZMANN EQUATION 


If the expression for W(z) is known, we can solve the 
mPB equation numerically to obtain the ionic density 
profiles. 



To test the simulation method developed above, 
we compare the results with a modified Poisson- 
Boltzmann (mPB) equation. The mPB equation is con¬ 
structed to account approximately for the ion-image and 
charge-charge correlations near an interface and is ex¬ 
pected to work well in the weak coupling limit. It was 
tested against MC simulation for a single dielectric in¬ 
terface with 1:1 electrolyte and was found to be very 
accurate. Therefore, we expect that for a slab geometry, 
a suitably modified mPB equation will also remain very 
accurate, allowing us to test the new simulation method. 

The electrostatic potential between two negatively 
charged dielectric surfaces satisfies the exact Poisson 
equation 

V^(j){z) = - [Qp+[z)-Qp-{z)] , (29) 

where is the mean electrostatic potential at a dis¬ 
tance z from the first plate (for simplicity it is placed at 
z = 0), p+{z) and P-{z) are the concentrations of cations 
and anions derived from salt and surface dissociation. In 
equilibrium, all the ions will be distributed in accordance 
with the Boltzmann distribution, exp(—, where Wg 
is the potential of mean force of ion of type s. We will ap¬ 
proximate Wg by the mean-electrostatic potential plus the 
correlation contribution, W{z). Suppose that there are 
Ab|_ cations and counterions, and N- anions per square 
Angstrom. Their distributions are then given by, 


FIG. 2. Representation of an electrolyte in the region < 
z < L — Vc used to calculate W{z). 

The correlational and charge-image contribution W{z) 
can be calculated approximately as follows. Consider 
ions (per A^), with charge Q and N- ions (per A^) with 
charge —Q, all with hydration radii rc, confined between 
two neutral dielectric walls with separation L, see Fig.[2l 
Due to the hardcore repulsion of ions from the surfaces, 
they are restricted to the region z > Tc and z < L — Tc- 
To keep the charge neutrality of this region, we intro¬ 
duce a uniform neutralizing background charge density 
—QNc/{L — 2rc), where A), = Ab|_ — N_. In the exterior 
regions z < 0 and z > L the dielectric constant is Co, 
while in the interior region it is e^,- The function W{z) 
then corresponds, approximately, to the energy penalty 
that an ion located at a distance z from one of the sur¬ 
faces feels due to asymmetry of its ionic atmosphere and 
due to its hard core repulsion from the wall. To obtain 
this potential, we calculate the Green’s function for a 
system of differential equations: Laplace equation, 

V2(/.(s, z) = 0 , (32) 

in the region with no electrolyte and a linearized Poisson- 
Boltzmann (LPB) equation, 

V^^(s, z) = K^(/)(s,z) , (33) 


p+(z) = 




f “ dz e-/5'3-/i(D-/3W(z) 
JVr 


PQ4>{z)-l3W{z) 

P-iz) = N-—j - 

Tr- 


(30) 


(31) 


in the region accessible to ions^, where (j){s, z) is the po¬ 
tential at position (s, z) in cylindrical coordinates and 
K = \/47rAB(A^+ -I- N-)/L is the inverse Debye length. 
LPB equation is used to account for the electrostatic cor¬ 
relations between the ions^. For 1:1 electrolyte, lineariza¬ 
tion of the Poisson-Boltzmann equation is justified since 
the ionic interactions are weak. Because of the azimuthal 
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symmetry of the problem, it is convenient to work with 
the Fourier transform of the potential, (/)(k, z) defined in 
terms of^. 


1 

(j){s,z) = -^ I dk kJo{ks)(j)(k, z) . 

Jo 


(34) 


For dielectric interfaces between hydrocarbons and water, 
eo/ctu ^ 1, so that to leading order we can set Co = 0. 
Taking into account the continuity of the electrostatic 
potential and of the normal component of the displace¬ 
ment field at z = 0, z = Tc, z = ZQ, z = L — rc and z = L, 
the Fourier transform of the electrostatic potential in the 
region 3, see Fig. [U can be calculated to be: 


MKz) = A3eP^ + B3e-P^ , 


(35) 


where k is the wave vector and p = \/kJ -\- k"^. The con¬ 
stants are given by: 




[(p -b + {p-k) 


_ p'\^2rcp+rck'} 


e-P^Q [{p + k)e^Lp+rck + (p _ fc)g2Lp-r,fej ^ ^ 

B3 = , 

/i =pcosh(fcrc) + /csinh(fcrc) , 

/2 =pcosh(fcrc) — fcsinh(fcrc) , 
/3=p(e2^P/2_e4’-=P/|) . 

Considering s —0 and z ^ zq, we find (now ommiting 
subscript 3): 

= 1^ r dkk{ [(p + A:)e^'-“P-’''='= + 

Jo J3 

{p - + lle-^P^Q+^^-P Up + k)e^LP+''->^ + 

h 

{p — fc)e^^P“’’‘'^]) -b June . 

(36) 


where June does not depend on the ion-plate distance zq 
and can be ignored. The work necessary to insert an ion 
at position z = zq from the interface, can be calculated 
using the Giintelberg^ charging process: 


Wi^Q) = ^ dkk{ 

Jo 

(p - k)e^^-P+^-'^] +^e- 
J3 


ke^p^Q [(p + fc)e2"-=P-’'-='= + 

J3 

2pZQ+2r,p ^■^^2Lp+r,k 


(p - k)e^^P-^‘^]) . 

(37) 




FIG. 3. Ionic density profiles for the cases: to = f-w = 80, left, 
and £0 = 2 and Cu, = 80, right. The plate charge density is 
—0.02 C/vaJ, the distance between plates is L = 50 A and the 
monovalent salt concentration is 10 mM. Symbols represent 
the simulation data and lines represent the solution of the 
mPB equation, Eq. (I29F Solid lines and circles represent 
cations, while dashed lines and squares represent the anions. 


profiles between two interfaces for the cases: (1) where 
£0 = 2 and Cw = 80, for which W{z) is approximately 
given by Eq. ( 133 , and (2) where eo = Sw = 80, for which 
W(z) is zero. The agreement between theory and simu¬ 
lations is very good, validating both methods. 



The interaction potential W{z) can now be used in the 

mPB equation, Eq. ([221), to account for the charge-image FIG. 4. Gouterion density profile for the cases aj = 3, left, 
interaction and the polarization of the ionic atmosphere. and aj = 5, right. The plate charge density is —0.02 G/m^, 

the distance between plates is L = 20 A. 


V. RESULTS 

We first show the comparison between the ionic con¬ 
centrations obtained using mPB Eq. (1^ and the results 
of MC simulations. In Fig. [3] we show the ionic density 


We next explore the effects of the dielectric hetero¬ 
geneity for strongly correlated systems. The valency of 
the counterions is modified to aj = 3 and aj = 5. The 
other parameters are kept the same as before. If ew = £ 0 , 
most of the multivalent ions adsorb to the charged wall 
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forming a strongly correlated quasi 2D one component 
plasma, see Fig. 01 However, in the case of large dielec¬ 
tric contrast between solvent and the confining surface, 
the counterions experience a strong charge-image repul¬ 
sion from the surface, see Fig. 01 This can significantly 
affect the interaction between charged dielectric bodies 
inside an electrolyte solution^!. Finally, in Fig. 01 we 
show the ionic distribution for 3:1 electrolyte between 
neutral dielectric surfaces. 



FIG. 5. Ionic density profiles for Co = 2 and 3:1 salt confined 
between neutral surfaces, for various salt concentrations. On 
the left, 100 mM of salt, while on the right, 500 mM. The dis¬ 
tance between plates is L = 50 A. Circles represent trivalent 
cations and squares monovalent anions. 


VI. CONCLUSIONS 

In this paper we developed a new simulation approach 
to study electrolytes and ionic liquids in a dielectric slab 
geometry. The method is easy to implement and is com¬ 
parable in time consumption with the regular 3D Ewald 
summation method. In the weak coupling limit we also 
presented a mean-field theory which allows us to pre¬ 
dict the ionic distribution between the dielectric charged 
plates. The agreement between both approaches, the¬ 
oretical and simulational, is very good, validating both 
methods. Examples of ionic density profiles for strongly 
correlated systems are also presented. Einally, the sim¬ 
ulation method developed here can be used to explore 
the interactions between colloidal particles with strong 
dielectric contrast. This will be the subject of the future 
work. 
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